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ABSTRACT: Perceptual studies suggest that the visual system uses the "rigidity" as- 
sumption to recover three dimensional structure from motion. Tillman (J984) recently 
proposed a computational scheme, the incremental rigidity scheme, which uses the rigid- 
ity assumption to recover the structure of rigid and nonrigid objects in motion. The 
scheme assumes the input to be discrete positions of elements in motion, under ortho- 
graphic projection. We present formulations of Tillman's method that use velocity infor- 
mation and perspective projection in the recovery of structure. Theoretical and computer 
analyses show that the velocity based formulations provide a rough estimate of structure 
quickly, but are not robust over an extended time period. The stable long term recovery 
of structure requires disparate views of moving objects. Our analysis raises interesting 
questions regarding the recovery of structure from motion in the human visual system. 
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1. Introduction 

An important source of three dimensional information is provided by the relative mo- 
tions of elements in the changing two dimensional image. The human visual system 
is capable of recovering structure from motion, under both orthographic and perspec- 
tive projection, and in the absence of all other cues to 3 D structure (see, for example, 
Miles, 1931; Wallach and O'Conneil, 1053; Bratmstein, 197G; Johansson, 1978; Tillman, 
1979). In studying the computation of structure from motion, one immediately faces 
the problem that the recovery of structure is underconst rained; there arc infinitely many 
3 D structures consistent with a given pattern of motion in the changing 2 D image. 
Additional constraint is required to establish a unique interpretation. 

Early perceptual studies suggested that the rigidity of objects may play a key role in 
the recovery of structure from motion (Wallach and O'Connell, 1953; Cibson and Cibson, 
1957; Green, 19G1; Johansson, 1975, 1977). Computational studies later established that 
rigidity is a sufficiently powerful constraint to derive a unique interpretation of structure, 
under a variety of viewing conditions. For example, Ullman and Kremlin (Ullman, 1979) 
showed that under orthographic projection, three views of four non coplanar points are 
sufficient to guarantee a unique 3 D interpretation (up to an unavoidable reflection about 
the image plane). In the case of perspective projection, Longuet Iliggins and Prazdny 
(1981) proved that the instantaneous velocity field and its first and second spatial deriva- 
tives at a point admit at most three different 3 D interpretations. Tsai and Huang (1981) 
showed that, with the exception of a few special configurations, two perspective views 
of seven points in motion are sufficient to guarantee a unique 3D interpretation. Wax- 
man and Ullman (1984) also addressed the uniqueness of the recovery of structure under 
perspective projection, basing their results on a kinematic analysis of continuous image 
flows. Additional theoretical results have been obtained for various classes of restricted 
motion, such as planar surfaces in motion (Hay, 19GG; Longuet Iliggins, 1984; Waxtnan 
and Ullman, 1985; Ulbuan, 1985; Negahdaripour and Horn, 1985), pure translatory mo- 
tion (Olocksin, 1980), planar or fixed axis rotation (Hoffman and Fliuchbaugh, 1982; 
Webb and Aggarwal, 1981; Bobick, 1983; Bennett and Hoffman, 1984a,b; Sugie and Ina- 
gaki, 1984), and translation perpendicular to the rotation axis (Longuet Iliggins, 1983). 
A review of the theoretical results regarding the recovery of structure from motion can 
be found in Ullman (1983). 

Prom theoretical studies of the structure from motion problem, it can be concluded 
that by exploiting a rigidity constraint. 3 D structure can be recovered from motion 
alone, using image information that is integrated over a small extent, in space and in 
time. These theoretical studies have also given rise to algorithms for deriving the rigid 
3D structure of moving objects (for example, Ullman, 1979; Longuet Iliggins, 1981; Tsai 
and Huang, 1981). Experimentation with these algorithms has revealed two important 
limitations. First, although it is possible in theory to recover structure from motion 
information that is integrated over a small extent in space and time, such a strategy 
may not be robust, in practice. A small amount of error in the image measurements 
can lead to very different solutions (Ullman, 1983). Second, most previous algorithms 
derive a three dimensional structure only when a rigid interpretation is possible, and 



otherwise do not yield any interpretation of si ruct tire or yield a solution that is incorrect 
or unstable. 

The lirsl observation above suggests that a robust algorithm for recovering structure 
should use motion information that is more extended in space or time. This conclusion 
is supported in recent computational studies by Nagahdaripour and Horn (1985) and 
Ullman (1984). Negahdaripour and Horn addressed the recovery of the motion of an 
observer relative 1 to a stationary planar surface. It was shown thai a robust recovery 
of both the observer motion and the orientation of the plane is possible when dense 
measurements of the spatial and temporal derivatives of image brightness are integrated 
over a large region of the changing image. Thus, consideration of motion information 
that, is more extended in space can lead to a stable recovery of structure. Bruss and 
Horn (I98,'>) also proposed an algorithm for recovering the motion of an observer relative 
to a stationary scene, which integral es motion information over an extended region of 
the image. The study by Ullman (1984), which will be developed further in this paper, 
demonstrated that a robust recovery of structure is also possible when motion information 
is integrated over an extended period of time. The extension in time can be achieved, 
for example, by considering a large number of discrete frames or by observing continuous 
motion over a significant temporal extent. 

With regard to the hnman visual system, the dependence of perceived structure on 
the spatial and temporal extent, of the viewed motion has not yet been studied systemati- 
cally, but the following informal observations have been made. Regarding spatial extent, 
two or three points undergoing relative motion are sufficient to elicit a perception of 3D 
structure (Borjesson and von Hofsten, 1973; Johansson, 1975), although theoretically the 
recovery of structure is less constrained for two points in motion, and perceptually the 
sensation of structure is weaker. An increase in the number of moving elements in view 
appears to have little affect on the quality of perceived structure (for example, Petersik, 
1980). Regarding the temporal extent of viewed motion, Johansson (1975) showed that 
a brief observation of patterns of moving lights generated by human figures moving in 
the dark (commonly referred to as biological mot ion displays) can lead to a perception of 
the 3 D motion and structure of the figures. Other perceptual studies indicate that the 
human visual system requires an extended time period to reach an accurate perception of 
3D structure (Wallach and O'Oonnell, 1953; White and Mucscr, i960; Green, 1961). A 
brief observation of a moving pattern sometimes yields an impression of structure that is 
"flatter" than the true structure of the moving object. (Tillman, 1984). Thus, the human 
visual system is capable 1 of deriving some sense of structure from motion information that 
is integrated over a small extent in space and time. An accurate perception of structure 
may, however, require a more extended viewing period. 

It was noted earlier that most algorithms for recovering structure from motion are 
unable to interpret nonrigid motions. There are, however, some exceptions to this. Ben- 
nett and Hoffman (1981b) studied the minimum amount of motion information required 
to derive a unique interpretation of the structure of a set of discrete elements undergoing 
nonrigid motion, when it is assumed that the elements are rotating about a fixed axis 
in space. Hoffman ami Flinchbaugh (1982) proposed an algorithm for interpreting the 
3 D motion and structure in biological motion displays. This algorithm decomposes the 



overall nonrigid motion inlo pairs of points that, arc rigidly linked and rotating in a plane, 
and triplets of points forming two hinged links that rotate in the same plane. Koenderink 
and Van Doom (1981; Koenderink, 1981) examined (lie class of bending deformations, 
which satisfy the physical constrain! that distances along the surface of the object, are 
preserved by the transformation. This class of deformations excludes any stretching or 
compressing of the object surface. In its current formulation, the method proposed by 
Koenderink and Van Doom for recovering the structure of bending surfaces requires that 
the surfaces ho complete, in contrast with other algorithms that are able to interpret the 
structure of isolated points in motion. To conclude, the algorithms discussed above for 
recovering the structure of nonrigid objects in motion all address restricted classes of 
these motions, such as lixed axis motion, planar motion, and bending deformations. 

The mechanism for recovering structure from motion in the human visual system 
appears not to be based strictly on the rigidity assumption. It is an everyday experience 
to perceive the structure and motion of deforming objects such as a flowing river, an 
expanding balloon, or a dancing ballerina. Perceptual studies reveal that, the human 
visual system can derive some sense of structure for a broad range of nonrigid motions, 
including stretching, bending and even more complex types of deformations (Johansson, 
1964, 1978; Jansson and Johansson, 1973; Todd, 1982, 1984). It is also the case that 
displays of rigid objects in motion sometimes give rise to the perception of somewhat, 
distorting objects (Wallach, Woisz and Adams, 1956; White and Mueser, 19G0; Green, 
1961; Braunstein, 19G2; Sperling et al, 1983; Uildreth, 1984a; Adclson, 1985). 

Iu this paper, wo focus on the recent work of Ullman (1984), which provides a 
more flexible method for deriving the structure of rigid and nonrigid objects in motion, 
and provides a natural moans for integrating motion information over an extended time 
period. This method makes use of the rigidity assumption, but in a different way from 
previous studies. The algorithm, called the incremental rigidity scheme, maintains an 
internal model of the structure of a moving object, which is continually updated as now 
positions of image elements are considered. The initial mode 1 ! may be flat, if no other cues 
to 3 D structure are present, or it may be determined by other cues available, for example, 
from binocular stereo, shading, texture, and perspective (Marr, 1982; Ballard and Brown, 
1982; Horn, 1985). As each new view of the moving object appears, the algorithm 
computes a new set of 3D coordinates for points on the object, which maximizes the 
rigidity in the tranfonnation from the current model to the new positions. In particular, 
the algorithm minimizes the change in the 3D distances between points in the model. 
The formulation presented by Ullman assumes the input to the recovery process to consist 
of a sequence of discrete frames, each containing a set; of discrete feature points. Through 
the process of repeatedly considering a now frame in the sequence and updating the 
current model of the structure of the moving features, the incremental rigidity scheme 
builds up and maintains a 3D model, and can be applied to both rigid and nonrigid 
objects in motion. Further details of the incremental rigidity scheme are presented in 
section 2 and in Appendix A. 

The incremental rigidity scheme has a number of advantages, from a computational 
perspective (Ullman, 1984): (1) because it integrates information over an extended time 
period, it provides a stable recovery of structure, particularly in the presence of error in 



the image measurements. (2) il allows deviations from rigidity, while always maintaining 
some model of Hie 3 I) struct lire of Hie object,, (3) it provides a natural means for 
interactions with other sources of 3 1) information, and (<1) empirical studies suggest, 
that the algorithm is able to recover the correct. 3 I) structure. The behavior of this 
algorithm is also consistent in several ways with human perceptual behavior (Oilman, 
198-1). 

In this paper, we develop Tillman's work further, in several directions. First,, in 
section 2, we present a continuous formulation of the incremental rigidity scheme, which 
uses velocity information at discrete points as input to the recovery process. In section 
3, we then examine in more detail, the behavior of the incremental rigidity scheme when 
presented with rigid objects undergoing rotation about a single axis in space, under 
orthographic projection. In particular, through computer simulations and a theoretical 
analysis, we examine the behavior of the discrete formulation as a function of the angular 
displacement between frames, and compare this behavior with that of the continuous 
formulation. Finally, in section 1, we present, both discrete and continuous formulations 
of the incremental rigidity scheme that, use perspective projection. Through computer 
simulations, we begin to examine the behavior of the perspective formulations, when 
presented with rigid objects undergoing both pure rotation about a single axis, and pure 
translation through space. 

The main conclusions of the paper are the following. The direct use of velocity 
information as input to the incremental rigidity scheme can provide a rough estimate of 
the structure of a moving object over a short viewing period, but is not sufficiently pow- 
erful to allow a detailed and robust recovery of structure over an extended time period. 
The computation of a stable long term solution appears to require the use of views of 
a moving object that differ significantly. This implies the need for a recovery process 
with memory of tin 1 past views, but this memory need not be extended indefinitely and 
continuously into the past. A small number of discrete views of a moving object are suf- 
ficient for recovering 3 D structure;. In the case of the incremental rigidity scheme, the 
use at every instant of a current, model of the 3 D structure of the object, and a present 
view that is sufficiently different from preceding views, can provide a robust recovery of 
structure. In the case of rotation of a rigid object about a single axis in space, both the 
rate of convergence of the algorithm to the final solution and the quality of the solution 
decrease as smaller angular displacements between viewed frames are considered. In the 
limit of the continuous formulation, the solution is no longer stable. The behavior of the 
perspective formulation of the incremental rigidity scheme is more complex than that 
of the orthographic formulation. We found that if the absolute' position of an object 
in space is known throughout the motion of the object, then the perspective formula- 
tion performs well, similar to the orthographic formulation. The results of computer 
simulations revealed a degradation in performance with smaller angular and spatial dis- 
placements between frames, but this degradation was somewhat more severe than in the 
orthographic formulation. Again, in the limit of the continuous formulation, the solution 
is no longer stable. Our analysis raises important questions* regarding the quantitative 
ability with which the human visual system can recover structure from motion; these 
questions are discussed in section 5. 



2. Discrete an<l Continuous Formulations of the Incremental Rigidity 
Scheme 

In this section, we first describe Tillman's (1981) original formulation oft lie incremental 
rigidity scheme, which assumes the visual input to consist of a. sequence of frames, each 
containing a number of discrete points that may correspond to identifiable features in the 
changing image. We then present a formulation thai, uses velocity information at discrete 
points in a continuously changing image as input, to the recovery process. The analysis 
in this section assumes orthographic projection of the scene onto the image plane. 

The motivations for considering a continuous formulation are threefold. First, on the 
basis of the results of computer simulations, llllman (1981) noted that, when analyzing 
objects undergoing rigid rotation, the convergence of the incremental rigidity scheme to 
the correct, solution was slower when smaller angular separations between frames wen- 
used. This suggests that (he scheme may perform bettor when successive views of the 
object differ significantly. We considered the limit, of arbitrarily close frames, both as 
a means of studying this phenomenon, and to determine whether a robust recovery of 
structure is still possible under these conditions. A second motivation is that recent 
work on the computation, of an instantaneous 2 D velocity field from the changing im- 
age suggests that a unique velocity field can be obtained for genera] classes of motion, 
exploiting a constraint on the smoothness of the velocity field (Horn and Schunck, 198.1; 
Hildreth, 1984a,b; Nagel, 1981). Ultimately, it may be useful to integrate the results of 
such velocity field computations with the recovery of structure from motion. A third 
motivation is that Ullman's formulation of the incremental rigidity scheme leads to the 
solution of a set of nonlinear equations. It is shown in Appendix A that the continuous 
formulation presented here leads to the solution of a set of linear equations. This makes 
a, theoretical analysis of the solution more accessible, and could in principle result, in a 
more efficient computer implementation. 

2.1 Ullman's Discrete Formulation 

The incremental rigidity scheme maintains and updates an internal model M(t) of the 
viewed object, which consists of a set of 3-D coordinates: M[t) = (xi(t), yi(t), Zi(t))- 
In this section, we assume orthographic projection (the case of perspective projection is 
addressed in section 4) onto the X — Y image plane, so that (xi(t),yi(t)) are the image 
coordinates of the i th point, and Zi{t) is the current estimate" of the depth at the i- -til 
point. (We assume a left handed coordinate system, with the positive z axis pointing 
away from the observer.) In Ullman's formulation, when no other 3D cues are present, 
the initial model M{t) at t = is taken to be flat; that is, ^(0) = (or some other 
constant value) for ? — !,... , n, where n is the number of points in motion. In principle, 
other initial configurations could also be considered. The theoretical analysis of section 3 
examines the long term stability of the incremental rigidity scheme, independent of the 
initial model of the structure of the moving points. 

Given a current model M(t) at time t, and the image of the moving points in a 
new frame at a later time t' , the problem is to compute a new model M(t') such that 



the transformation from A/(/) to M(f') is as rigid as possible Since x L (l') and ?/,(/') are 
known, this requires the computation of Hit- unknown depth values ~ t (/'). (It is assumed 
that the correspondence between points in the two successive frames is known.) The new 
depth values are computed as follows. Let /,;_/(/) denote the distance between points i 
and j at time t. To make the transformation as rigid as possible, the values Zi(t') for 
the new model are chosen so as l.o make /;_/(/) and hj(f') as similar as possible. For this 
purpose, llllman defined a measure of the difference between /,;/(<) and /,/(/') as: 



(UAt) ~ hA?)) 
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dikjitYkjit')) = ^1-L^-LL-, (2.1) 

and formulated the recovery of structure as (lie computation of Zi(f') that minimize the 
following' overall deviation from rigidity: 

D(tj') = Td(l ij (l)J ij (t')). (2.2) 

After tlie values Zi(t') have been determinod using this minimization process, the new 
model M(t') — [xi[t') ,'!Ji{t') , Zi(t')) becomes the current, model. A new frame is then 
registered and the process repeats itself. In this way, the scheme maintains rigidity by 
keeping the total distances between points in the model as constant as possible. The 
motivation for the cubic factor in the denominator of Eq. (2.1) is that; the nearest 
neighbors to a given point are more likely to belong to the same object than distant 
neighbors, so that a point is more likely to move rigidly with its nearest neighbors. The 
if jit) factor diminishes the influence of distant points in the recovery process. 

It should be noted that in the case of orthographic projection, only relative depth 
values, Z{(t) — Zj{t). can be recovered, rather than absolute depth values, because under 
this form of projection, the image of a given object does not change with its absolute 
depth. In addition, 3-D structure is determined only up to a reflection about the image 
plane, since the orthographic projection of a rotating object, and its mirror image rotating 
in the opposite direction, coincide!. 

2.2 The Continuous Formulation 

A continuous formulation, which uses velocity information at discrete feature points, 
can be developed as follows. Assume again that there always exists an internal model 
M[l) — (xi(t),yi(t),Zi(t)). Assume also that the image velocities ii(t) and iji{t) are 
known. The problem is then formulated as the computation of the z components of 
velocity, z.i{t), that minimize the total continuous change in the distances between the 
points. The general form of the measure of overall deviation from rigidity is given by: 

*M0 = £ 'MM'))- (2-3) 

In our.analysis, we consider different possibilities for the measure, d c (hj {(■))• The theo- 
retical development of section 3.2, for example, considers the behavior of the incremental 



rigidity scheme, as the frames become inliiiilesimally close, using the following discrete 
measure of the change in I. lie distance between pairs of points: 

rf(/, V (0^.v('')) = (^-(0 " JiyC')) 2 - (2.4) 

The derivation of the continuous measure of rigidity is outlined in section 3.2. 1. It, results 
in the following expression for D r (t), as a function of the coordinates and velocities of 
the points (the arguments, /, have been omitted for simplicity): 

D e {t) = Y, ^ -' X ^ ii - *y) + (»' ~ yj)(in -■ w) + (- - *;)(* - ^)) 2 - ( 2 - 5 ) 

This particular measure of rigidity was used in the theoretical study for analytic sim- 
plicity. The computer simulations of the continuous formulation use the measure of 
change in the distance between points given again by the limit of (l,j(t) - lij(t')) 2 , for 
infinitesimally close frames, which is: 

<MM0) = ta(*)) 2 . (2-6) 

In terms of the coordinates and velocities of the points, this yields the following overall 
measure of deviation from rigidity: 

cW " ^ (*.- - xy) 3 + (:</; - yy) a + (*» - ^) 2 " l ' J 

(Eqs. (2.5) and (2.7) use slightly different measures of rigidity, but, serve the same role 
as measures of overall changes of rigidity for the continuous formulation. In the present, 
paper we do not, use different notations for these measures, as it will be clear from the 
context which measure is used. More specifically, Eq. (2.5) is used in the theoretical 
analysis of section 3.2, and Eq. (2.7) is used in the computer simulations of section 3.1.) 
In other respects, the continuous formulation is similar to the discrete formulation. A 
model of the structure of the moving points is built up by continually taking into account 
new velocity information over an extended time period. Again, because orthographic 
projection is used, only relative velocities, i,({) - z 3 -(t), can be recovered. This can 
clearly be seen in Eqs. (2.5) and (2.7), in which the coordinates of the points and their 
time derivatives all appear in differences between pairs. Further details of the continuous 
formulation are presented in section 3 and Appendix A. 

The analyses presented in this paper mainly consider single rigid objects in motion, 
which are compact in the sense that the internal distances between pairs of points do not 
differ much from one another. In tins case, the additional /?• factor of Eq. (2.1) has little 
influence on the behavior of the algorithm, so we omitted it in our theoretical analysis and 
computer simulations for the sake of simplicity. In general, however, a proper weighting 
(and not, necessarily a cubic factor) of the influence of different distances among points 
is necessary for a better performance of the algorithm. 
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3. Positions vs. Velocities as Input to the Recovery of Structure 

We staled earlier thai, on the basis of computer simulations, Ullmaii (lOS-l) observed 
that when analyzing objects undergoing rigid rotation, the convergence of the incremen- 
tal rigidity scheme to the correct solution was slower when smaller angular displacements 
between frames were used. In this section, we analyze this phenomenon from a theoret- 
ical perspective, focusing on the long term stability of the computed 3 D model. We 
first, examine the behavior of the continuous formula! ion, which uses a current model 
and measured image velocities at discrete points as input to the recovery process. We 
then turn to the discrete formulation, which uses the discrete positions of discrete points 
as input, and examine the long term stability of its solution as a function of the angular 
displacement between frames. Our main conclusions are the following. First, for the 
particular class of motions considered, the discrete formulation always yields a 3 D so- 
lution that converges asymptotically to the correct solution, but the rate of convergence 
varies with the angular displacement. The rate of convergence increases with increasing 
angular displacement up to a maximum, ami then decreases with further increases in 
this displacement. The position of this maximum depends on such factors as the type of 
motion and geometric structure of the points. 

Although the orthographic projection is in general not physically valid, it, is used here 
because it, allows a simpler formulation of the problem, and is therefore better suited to 
theoretical analysis and computer implementation. It allowed us to gain a deeper insight 
into the nature of the phenomena studied. We nevertheless implemented the equivalent 
perspective formulation and confirmed that, the basic results remain valid under this 
formulation. The use; of perspective projection enables the recovery of the structure 
of objects undergoing pure translation, which was not possible under the orthographic 
projection. In the case of translation, the rate of convergence of the computed 3 D model 
to the true structure also increases with increasing spatial displacements between frames. 

Before presenting the results of the theoretical analysis, we illustrate the behavior 
of the discrete and continuous formulations through the results of computer simulations. 
We show that the continuous formulation yields an initial fast, convergence to a close 
approximation of the true structure of the moving points, but then oscillates over a large 
range. It consequently does not yield a stable long - term recovery of structure. 

3.1 Observations from Computer Simulations 

In this section, we briefly illustrate; the behavior of the discrete and continuous formula- 
tions of the incremental rigidity scheme, for the special case of rigid rotation of a small 
set of discrete points about the vertical axis. (For details of the computer implementa- 
tion, see Appendix A.) In the case of the discrete formulation, we examine the rate of 
convergence of the algorithm and the quality of the solution that it yields, as a function 
of the angular displacement between frames. We then compare its performance with that 
of the continuous formulation. In all of the examples presented here, the input, consisted 
of a set of five points in space. The first point is assumed to lie at the origin of a co- 
ordinate system that is displaced from the viewer along the line of sight. The position 
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of this first point, is constant, throughout the motion of the points. Tin 1 coordinates of 
the remaining four points were chosen randomly. Fig. 1 illustrates a typical set. of five 
points, showing their projections onto the X Y plane (Fig- la) and the X Z plane 
(Fig. Ih). In the simulations, the projected positions and velocities of the points were 
computed analytically, rather than measured from real image sequences. 
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Figure 1. A set of five points with random coordinates is projected onto (a) the X — Y plane, 
and (b) the X — Z plane. 



Fig. 2 illustrates the behavior of the discrete formulation of the incremental rigidity 
scheme, as a function of the angular displacement between frames. Each figure shows 
a birds' eye view of the set of rotating points (that is, their projection onto the X — Z 
plane), with filled circles representing the true positions of the points and open circles 
representing the structure computed by the algorithm. Fig. 2a shows the true positions 
of the points and the initial model at time t = 0. The initial model is assumed to be 
flat. Figs. 2b and 2c show the true and computed configurations of points after 120° of 
rotation. This final position was reached by taking three steps of 40° (Fig. 2b), and 12 
steps of 10° (Fig. 2c). It can be seen that the use of a smaller number of more disparate 
views yields a 3D model that is closer to the true solution. As noted in Appendix A, a 
steepest descent minimization algorithm was used for most of our computer simulations. 
We also analyzed a small number of examples using an exhaustive search algorithm to 
find the minimum solution, and again found the accuracy of the solution to vary with 
the angular displacement, between frames. We therefore believe this to be a fundamental 
behavior of the incremental rigidity scheme that is not simply a consequence of the 
particular algorithm used to implement the scheme. This observation is supported further 
by the theoretical analysis of the next section. 

In Fig. 3, we show a series of graphs that illustrate 1 in a different way, the behavior 
of the discrete formulation of the incremental rigidity scheme, as a function of angular 
displacement. In this case, the set of points shown in Figs. I and 2 were rotated by 
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Figure 2. (a) The true configuration of five points (filled circles) is compared with the initial 
configuration of the points in the model (open circles) at time t — 0. The projection is onto the 
X — Z plane, (b) The comparison between the true and computed positions of the points after 
three stops of 40°. (c) The comparison between the true and computed positions of the points 
after 12 steps of 10°. 



four full revolutions. Each of the graphs show the error between the true and computed 
structures, as a function of time. In particular, the following quantity is plotted: 



y,(<hj - lij-y- 



*'tf 



(3.1) 



where d{j is the correct 3 D distance between points i and j in the object, and lij is 
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Figure 3. Graphs of the error. in the internal distances between points in the computed 3D 
model, as a function of time. The points are rotated 4 full revolutions, in steps of (a) 40", (b) 
20", (c) 10°, (d) 5° and (e) 1°. (f) The graphs in 3a 3e are shown superimposed. 
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the corresponding distance between points *' and j in (lie computed model. Tin? graphs 
shown in Figs. 3a through 3e correspond to rotations with angular displacements of 40", 
20", 10", 5" and 1", respectively. In Fig. 31', the five graphs are shown superimposed. 
Again, it can he seen that (lie rale of convergence and quality of the solution improves 
with larger angular displacements. For a particular total angular extent, the error in the 
computed model decreases with increasing angular displacement between frames. In the 
case of 40'' displacements, the algorithm converges asymptotically and monofonically 
to the final solution. For smaller displacements, the convergence is no longer strictly 
monotonic, but is still essentially asymptotic, toward the final solution. Fig. 4 shows 
the same set of graphs superimposed, but with the error plotted on a log scale. The 
convergence of the solution is now essentially linear, with varying slopes, suggesting that 
the actual convergence is exponential. 
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Figure 4. Graphs of the error in the internal distances between points in the computed 3D 
model, as a function of time, plotted on a log scale. The points are rotated 4 full revolutions, 
in steps of 40°, 20°, 10°, 5° and 1°. The graphs are shown superimposed, with the angular 
displacements indicated above each graph. 



Fig. 5 illustrates the behavior of the continuous formulation of the incremental 
rigidity scheme. The same set of points used previously was again rotated about the 
vertical axis and the 3D model was computed at infmitcsimally closely spaced times, 
using the instantaneous velocities projected onto the image (see Appendix A for details). 
In Fig. 5a, we compare the true positions of the points (filled circles) with the best 
solution (open circles) obtained over 10 full revolutions of the points. Although the model 
is quite close to the true structure at this position of the points, the solution oscillates 
significantly over an extended time period. A graph of the error in the computed model 
over the 10 revolutions is shown in Fig. 5b. There is an initial fast convergence toward 
the true structure of the points, but the algorithm then oscillates with high amplitude 
away from the true solution. When the error in the model is high, depth reversals often 
occur. Fig. 5c shows an example of the true and computed structures at a time of 
complete depth reversal. Such reversals were also observed with the discrete formulation 
of the scheme, although rarely. 

Fig. G illustrates the typical behavior of the discrete 1 formulation of the incremental 
rigidity scheme, averaged over 10 configurations of live points. For each of the conligura- 
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Figure 5. (a) The best solution (open circles) obtained over 10 revolutions of the points is 
compared with the true positions (Oiled circles), for the case of the continuous formulation, (b) 
The error in the internal 3D distances between points in the model, as a function of time, (c) 
A complete depth reversal between the true and computed structures. 



tions, llio first point, was placed at the origin oft.hr coordinate system and the coordinates 
of the remaining four points were chosen randomly. Tin- set of points was then rotated 
by three discrete angular steps, with the size of the angular displacements ranging from 
1° to 90°. The initial ,'5 D model for the points was Hat, and a new model was computed 
for each of the three discrete positions of the points. After the third step, we computed 
the following measure of the absolute error in the internal distances between points: 



E 



* o 
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(3.2) 



whore r^j is the (rue 3 D distance between points i and ;/, /,y is the distance between 
points i and j in the computed model and n is the total number of pairs of points. 
This measure was chosen because it expresses an average of the error in the model 
relative to the true structure. (Note that a lower value for this measure corresponds to 
less error in the computed model.) We also considered other measures of error in the 
distances between points and in the actual depth values, and found the general behavior 
of the algorithm to be the same under different measures. The above error measure was 
averaged over the JO configurations of points, and is plotted in Fig. G, as a function of 
angular displacement. It. can iirst be seen thai in general, the error after three discrete 
stops of the algorithm varies with the size of (.he angular displacement between frames. 
There is a steady improvement in performance as the displacement increases, to about 
50°, followed by a dogrodation for increments of GO , and a steady decrease in performance 
from 70° to 90°. The degradation in performance for an angular displacement of G0° was 
common; 8 of the 10 configurations of points exhibited this behavior. Tins degradation 
may bo a consequence of the symmetry between the initial and final views, which are 
rotated 180° from one another. In general, the convergence of the algorithm for three 
discrete steps degrades significantly for smaller angular displacements. This result is not 
surprizing, in that there is very little change in the discrete views for such small angles 
and a reduced total angular extent. The deterioration for large angles probably occurs 
because at, 90 c , the number of views available to the structure from motion process is 
reduced to two. 




70° 80° 
Angular Displacement 



Figure G. Absolute error in the internal distances between points in the 3 D model computed 
by the incremental rigidity scheme. The graph shows the error after three discrete views of the 
rotating points are considered, averaged over 10 random configurations* of tive points. 
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Fig. 7 ill ust rales (he error in the computed model, as a ('miction of angular dis- 
placement, for I lie case in which the overall rotation of the points was kept essentially 
constant,. The same set, of 10 random configurations oT points was rotated by discrete 
angular steps, with the si/.** of the steps varying from 10" to 90''. In this case, each 
configuration was rotated by a total of (approximately) 180' and 300", and the same 
measure shown in Eq. (3.2) was then computed. This measure was again averaged over 
the 10 configurations of points, and is plotted as a function of angular displacement in 
Fig. 7. Fig. 7a shows the data for the case of 180" of rotation. Note that for the angles 
40°, 50", 70" and 80", there are two points plotted, corresponding to multiples of the 
angular displacement that are just less than and greater than 180". The graph that is 
superimposed on the points passes between the pairs of points for these angles. Fig. 7b 
shows the same data for the case of 300 " of rotation (for the case of angular displacements 
of 50 and 70", the points were rotated by a total of 350 ). Also shown in Figs. 7a and 
7b are (lit- average errors in the solution (hat wore derived by the continuous formulation 
of the algorithm after 180 and 300 of rotation of the 10 configurations of points. These 
two data points are indicated by the stars along the ordinate of the graph. The main 
observation to be made is that when the points are rotated by a constant, total amount, 
there is still a strong dependence of the rate of convergence of the solution on the size 
of the angular displacement, between frames. There is again an improvement, in conver- 
gence rate as this angle increases, up to about 50°, followed by a slight worsening for 
an angle of 60°, then improvement again for 70", followed by a steady worsening to 90°. 
The degradation in performance for an angular displacement of 00° was again common, 
occuring for 7 of the 10 configurations of points. The deterioration of the convergence 
rate with decreasing angular displacements is not obvious, because in spite of the fact 
that the changes between consecutive frames arc smaller, there are many more frames 
altogether. In addition, while the continuous formulation provides a good estimate of the 
true solution after 180° (see Fig. 7a), the solution then degrades significantly, providing 
a relatively poor solution after 300° of rotation (see Fig. 7b). 

To conclude, the results of computer experiments with the discrete formulation of 
the incremental rigidity scheme show a clear dependence of the behavior of the computed 
solution on the size of the angular displacements between frames. In the limiting case 
of the continuous formulation, there is an initial fast convergence of the solution toward 
the true solution, followed by a substantial oscillation of the solution. The next section 
presents a theoretical analysis that attempts to explain this phenomenon. 

3.2 Analytic Study of Convergence Properties 

In this section wo outline the main conclusions of our theoretical analysis. The purpose of 
this analysis is not a general study of the behavior of the incremental rigidity scheme, as 
such an analysis would be too cumbersome. Rather it concentrates on a formal analysis 
of the convergence properties of the algorithm, for a family of examples that we believe 
are relevant, to the general recovery of struct lire from motion. We emphasize the concepts 
raised by this analysis, and comment on their generality. This section considers a subset 
of those examples analyzed in the computer simulations of section 3.1. 
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Figure 7. Absolute error in the internal distances between points in the 3-D model computed 
by the incremental rigidity scheme. The graph shows the error after a total of (a) 180° of 
rotation, and (b) 300° of rotation. The errors are again averaged over 10 random configurations 
of five points. The stars along the ordinate indicate the data points for the average errors of 
the continuous formulation after 180° and 300°. 



To simplify the theoretical analysis, we consider in this section a measure of overall 
deviation from rigidity that is somewhat different from the one used by Tillman, and it 
is the following: 

D d (t,t') = £(/?.(*) - l^(t'))\ (3.3) 

i,3 
This new measure uses the difference between the squares of the distances, ijAi), rather 
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Mian the difference between (he /{/(/)'s themselves, thus avoiding I he use of square roots. 
This measure also does not contain Hie cubic factor, lf.;(t). in the denominator, which was 
included in the measure proposed by Ullman as a means of reducing the contributions 
of distant points relative to nearby ones. In this theoretical analysis, we mainly consider 
configurations of points that are compact, in the sense that the internal distances between 
pairs of points do not. differ much from one another. In this case, the cubic factor seems 
not to be important. 

The main problem that this theoretical analysis addresses is the long term stability 
of the solutions obtained by the incremental rigidity scheme, as a function of the angular 
separation between successive frames, for the case of rotation of rigid objects about ;i 
single axis in space. In the previous section, we observed the variation of convergence 
rate with angular displacement, in the results of computer simulations. In the present 
analysis, we examine this convergence phenomenon through a study of the stability of 
the incremental rigidity scheme. In particular, we analyze the behavior of the internal 
model under small perturbations, when it is near the true solution. We first examine the 
limiting case, where the frames are arbitrarily close to one another, that is when t' —* t. 
The analysis of this case is presented in sections 3.2.1 and 3.2.2. We then explore the 
convergence properties of the discrete formulation in sections 3.2.3 and 3.2.4. Finally, in 
section 3.2.5 we summarize the theoretical properties of the two schemes, as a function 
of angular displacement. 

3.2.1 The Continuous Formulation 

In the case of the continuous formulation, lij{t') can be determined to a good first order 
approximation by /,;,(£) and its derivative, lij{i), that is: 

l ij {t')^h j {t) + l l . J {t)x{t' -t). (3.4) 

As described earlier, when such an approximation is used, rather than computing the 
depth values Zi{t') that maximize rigidity, we can compute the z components of velocity 
Zi(t) that do the same task. It can be shown that in the limit, as t' — > t, using Eq. (3.3), 
and the definition of lij, the quantity to be minimized is: 

D C {1) = J2U Xi ~ x N Xi ~ *i) + (»«' ~ %'M& ~ Wj) + (*•' ~ z i)i k i - z j)f- ( 3 - 5 ) 

Thus, given the model at time t , and the x and y components of velocity in the image, 
we compute the temporal derivatives of the depth values that minimize D c (t) given by 
Eq. (3.5). This is the continuous formulation that we use in the theoretical analysis of 
the convergence properties of the incremental rigidity scheme. 

Note again, that the quantities that can be computed are not the absolute 1 values of 
z,;(t) but the relative values z t (t) — Zj(t). It follows that without loss of generality, we can 
set (i(),?/o, 20) at the origin of the coordinate system, that is (*o,'(/o, i'o) — (0,0,0). At 
every instant, the remaining n — 1 points are then given relative to this first point. Thus 
the number of independent variables in the model is n — 1. Introducing the following 
notational simplification: 
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Oij -" (x.i - *j-)(ii • ij) I (?/,; - yj){yi - ?//), (3.0) 

Eq. (.5.5) c;ui then he rewritten a.s follows: 

n 2 n 1 n 1 

'MO = E X>* f te - *;K*«- - ^-)) 2 + X><>y + My) 2 - ( 3 -7) 

» i y ■« y j 

As described above, we are looking for the Zi(t) that minimize D r (t). The necessary 
condition for such a minimization is that the partial derivatives of D,.(t) with respect to 
the Zi(t) are zero: 

§^ r , ,, 8) 

In section 3.5 we show that Eqs. (3.8) represent a set of n I linear equations with n — 1 
variables Zi(t), which generally have a unique solution. It follows that because D c {t) > 0, 
thus being bounded from below, the £.,•(/.) that satisfy Eqs. (3.8) also minimize D,.(t), so 
that these equations generally represent a sufficient condition for the minimization. 

The z.i(t) that satisfy Eqs. (3.8) are expressed in terms of Z{(t), thus representing 
a system of n - 1 differential ecjuations with n — 1 variables. This system, however, is 
generally difficult to solve explicitly, as it is nonlinear. The only solutions that we were 
able to verify by straight substitution are the true motion of the object and its depth 
reversal (which are equivalent under orthographic projection). In the next section wc 
present a stability analysis of the true rigid motion solution; that is, wc examine the 
stability of the algorithm when its solution is close to the true solution. (The computer 
simulations address the full convergence behavior of the algorithm.) We concentrate on 
rigid motion, because we feel that nonrigid transformations may disrupt the stability 
of the solutions, introducing instabilities due to peculiar types of motion. For example, 
the internal model may be too sluggish in adjusting itself when fast changes of structure 
occur. By considering rigid motion in this analysis, we address more fundamental aspects 
of the theory itself. 

3.2.2 Stability Analysis of the Continuous Formulation 

The idea of a stability analysis can be stated as follows. Suppose that at a given instant t^, 
the computed 3 D model is very close- to the true solution, that is ^(<o) = Zi(to) + fi{to) 
where Zi(t) is the true depth value at point i. Because the system is perturbed at to, we 
expect it in general to be perturbed for every t > i () , that is Z{(f.) — Z{(t) + t{(t). The 
system is said to be asymptotically stable if the following is true: 

liuif»(0 = 0, (3.9) 

for every i. If the ti{i) remain bounded, but do not converge to zero, the system is 
defined to be weakly stable. If, however, lim t >00 (;(/-) — oo for some i, the system is said 
to be unstable. 
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If we suppose thai /,(/,) is small enough for every / > /, ( ). Mien we can make the 
following first, order approximation of ICqs. (3.8) in terms of <,(/.) and <',;(/) (I lie arguments, 
/,, have been removed for simplicity): 



y! f D 2 D, ( , i) 2 D c . 

f-J \ifi l i)z/ J dz-Oij 1 - 



0. (.3.10) 



A derivation of Eqs. (3.10) is given in Appendix B. Eqs. (3.10) can be readily solved for 
(j(t) in terms of (j(t). For notalional simplicity, we write- this solution in matrix form: 



i) 



<) 2 Dr 



, « .,. -, Af -, (3.11) 

i)Zi()z.j\ [OziOzj] 

where i = (< i , ..., < „ [ ) r , \<) 2 D,./i)zii)zj\ is the matrix whose dement at row i and column 
7 is <) D c ji)zi<)zj and \i) D c /i)zii)zj\ is the inverse of the matrix wliose element at 
row i and column j is <)~ ' D c jdz.,()zj, provided it exists. The matrix A is defined to be 
~~[() 2 D c /dzi(')zj] l [i) 9 'D,./i)zi()zj}, and is evaluated at the true solution itself at every 
instant. 

Using the Eqs. (3.11) to study the stability of the system, it is not possible to prove 
that the system is unstable. To do so requires a proof that r is unbounded as t — » oo, 
but in this case the first order approximation used to derive Eqs. (3.10) no longer holds. 
It is still possible to determine, however, whether or not the system is asymptotically 
stable. Indeed, we will show in this section, that for many types of motion, the continuous 
formulation of the incremental rigidity scheme is not asymptotically stable. The results of 
computer simulations provide evidence that the formulation is, in general, not unstable, 
thus being weakly stable. 

The Eqs. (3.10) represent a system of ordinary linear differential equations, which 
are used in conjunction with Eq. (3.7) where D c (t) is defined. Note therefore, that A 
has time dependent components. How do the different types of motion determine A? A 
generalized motion of a rigid body can be described instantaneously as a rotation about 
a iixed axis in space plus a translation of the origin of the coordinates. We noted before 
that under orthographic projection, all that, can be determined are the relative depths of 
the points and not their absolute values. Also, in the same case, the only relevant data 
to the problem of recovering structure 1 from motion are the relative image coordinates 
(terms d{j in Eq. (3.7)). Translation does not change any relative distances, as it changes 
the positions of all points by the same amount. Thus oidy rotations need to be considered 
in this problem. 

As mentioned earlier, the matrix A in Eq. (3.11) is in general time dependent. Thus 
this system cannot be solved by the standard method of characteristic values, available to 
systems with constant coefficients. Furthermore, the general rotation can have variable 
angular velocity, and the axis of rotation can change over time, making Eq. (3.11) very 
difficult to integrate analytically. 

It is not necessary, however, to solve Eq. (3.11) in order to reach certain conclusions 
regarding the nouconvergence off to (). For example, let. l,(t) be n— 1 arbitrary solutions 
of Eq. (3.11) with initial conditions r •(()) specified. Then the n - 1 X n - 1 matrix E{t) - 
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[< I (/.), ...,( n ,(/.)] satisfies the Liouviile Jacobi formula. (Yakubovich and Slarzhinskii, 
1975): 

IM{E{t.)) -■- Ihl.(E(()))cxp ( i Tr(A(t'))dA . (3.12) 

It follows (hat a necessary condition for liin t >QO E(/.) -•- 0, that is, for Hie system to be 
asymptotically si able, is: 

y»do 

/ Tr{A{t'))<W = -oo. (3.13) 

We found that condition 3.13 does not. hold for many interesting types of motion. 
For example, for any three point configuration, rotating with arbitrary angular velocity 
about an axis perpendicular to the viewing axis, we can derive the following relationship, 
by using Eqs. (3.7) and 3. J I (details of the derivation are given in Appendix C): 

Z 2 - Z{Zi + z t 

where Z{, Z2,?i, Zi denote the z coordinate's and velocities evaluated at time t. This trace 
is periodic and has zero integral over the cycle, for rotations of fixed angular velocities 

w: 



■A) 



-2tt/w 

Tr{A{t))dt - 0, (3.15) 

'o 

a result that contradicts Eq. (3.12). Thus, for any three point configuration have this 

type of motion, the internal 3 D model computed by the continuous formulation will not 

converge to the true structure. 

The fact that this result is derived for three point structures is nontrivial. In fact, it 
is shown in the analysis of the discrete formulation, later in this section, that converging 
models can be constructed for three point configurations under certain conditions. A 
detailed analysis for structures with a higher number of points is in general very cum- 
bersome, but we verified that Eq. (3.15) holds for special four point configurations, 
such as those that when viewed along the rotation axis, appear as squares, rectangles or 
trapezoids. 

We also note that having the rotation axis perpendicular to the viewing axis is the 
best situation, as far as convergence is concerned, because depth motion is lost as the 
rotation axis is slanted towards the viewing axis (that is, away from the image plane). 
In particular, note that in the limiting case, whore the rotation axis is parallel to the 
viewing axis, no motion in depth occurs whatsoever, and so the structure cannot be 
recovered from relative motion. 

Is the instability of the continuous formulation due to properties of the projections of 
the object at particular positions in its revolution, which somehow convey less information 
to the recovery of structure from motion (for example, when many of the points belonging 
to an object have little motion in depth)? In order to check this hypothesis, we explored 
the convergence of the internal model when the object was oscillated back and forth 
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over a small angular extent, around different positions. II was found again, by checking 
the validity of Eq. (3.15), that regardless of position, the internal model is unable to 
converge to (lie true structure, thus giving a negative answer to the question raised. This 
suggested to us that the reason underlying the instability of the continuous formulation 
may be based on the smallness of the angular displacement, between consecutive frames, 
and not on the structure of the frames themselves. Pursuing this idea further, the next 
sections address the convergence behavior of the discrete (and more general) formulation 
of the incremental rigidity scheme, focusing on the effects of the size of the angular 
displacements between frames. 



3.2.3 The Discrete Formulation 

The continuous formulation was derived from the approximation of Eqs. (3.3), as the 
angular displacements between the consecutive frames used in the structure from motion 
process become small, and it was expressed in Eq. (3.5). If we do not make such <ui 
approximation we have a more general formulation, which is however discrete. By being 
more general, this formulation allows a study of the incremental rigidity scheme as a 
function of the angular displacement between frames, and enables us to verify whether 
the total instability found in the continuous formulation is due to an "artifact" of the 
approximation made, namely an infinitesimal angular displacement between frames, or 
rather is due to a continuous process in which the system becomes less stable with smaller 
angles and eventually becomes unstable when the displacement is small enough. 

In order to stress the discrete nature of the general formulation, we first rewrite Eq. 
(3.3), using as independent variable t = kr, where k is an integer and r a fixed interval: 



D d (kr) = 2 [(xi(kT) - .xy(fcr)) 2 + ( yi (kr) - yj (kr)) 2 + (zi(kr) - z^kr)) 2 
id 

- ( Xi ((k + l)r) - x,-l(k + l)r)) 2 - (y t -((* + l)r) - y 3 {{k + l)r)) 2 < 3 - 16 ) 

- (*,-((* + l)r) - *,■((* + l)r)) 2 ] 2 . 

Note that as in the continuous formulation, the quantities of interest for the theory are 
relative rather than absolute, i.e. (zi((k + l)r) — Zj((k + l)f)). Thus again, without loss 
of generality we can set (xo,yo, 2 ( j) = ((),(),()), and study the behavior of the other n — 1 
points relative 1 to it. For notational simplicity, let: 



b i3 (kr) =( Xi ((k -I- 1)t) - xj((k h l)r)) 2 + (».-((* + l)r) - Vj ({k + l)r)) 2 

- ( Xi (kr) - Xj(kT)) 2 - ( yi [kr) - yj (kr))\ l °' ' 

We now rewrite Eq. (3.16), using this notational simplification and without explicitly 
specifying r: 
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In the present, formulation wo arc looking for (lie Zi[k-\- I) that minimize D,{(k). As 
in the continuous case, the necessary and suflicioi.it condition for this minimization is: 






(3.19) 



which represents a set of n I nonlinear o<|iiations implicitly relating the Zi{k) and 
Zi{k-\ I). The only instance in which we have been able to derive a full analytic solution 
for these equations is in the two point, case. This case is a degenerate one, because for 
two points the 3 D structure is not determined uniquely by any number of views. It is 
important to discuss this case, however, as it provides a relevant insight into the theory, 
which is discussed further in section 3.2.5. Eqs. (3.18) and (3.19) in this case reduce to: 



[b Ui {k) + z\{k + 1) - zl(k)] Zy (k + 1) = 0. 



(3.20) 



The only non trivial solution for this equation occurs when the term inside the brackets is 
zero. But using the definition of bm, this condition essentially implies that the distances 
between the two points is constant. In other words, the distance 1 between the points in 
the initial internal model at k — will tend to remain the same throughout the entire 
motion even if it is wrong. The model will only correct itself if its initial guess for the 
distance between the points is small. In that case, the internal model will be forced 
to change because the distance between the points in the projected plane will become 
larger than the distance in the model. In this case, numerical calculations show that 
the model expands, and even has a small overshoot such that the distance between the 
points becomes a little larger than the true distance, and then stays at this condition 
forever. 



3.2.4 Stability Analysis for the Discrete Formulation 

The stability analysis for the discrete formulation begins with the same general argu- 
ment as for the continuous case. Suppose that at. k — 0, the model is very close to the 
true solution, i.e. Zi(0) = £;(0) H- 1 ,(()), where £';(&) is the true depth value. We exam- 
ine whether or not the perturbation at later times converges to zero, that is, whether 
lim/ c >oo<i{k) — 0. If this is the case, and the (,:(&) always remain small, then we can 
write an equation similar to Eq. (3.11) for the discrete formulation: 



(_{k + i) 



iPD d 



Dziik-^^dzjik + l) 



<) 2 n d 



dz.i{k+ l)dzj(k) 



c{k) = B{k) L (k). (3.21) 
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Again, by the same arguments given in the continuous cases, we cannot use Eq. (.'5. 21) to 
prove Uie system to he unstable, but only verify whether its convergence is asymptotic. 

In this analysis we concentrate on three point configurations, because as in the con- 
tinuous for.iiiulci.Uou, an analytic study with a larger number of points is too cumbersome. 
In flie results of computer simulations, however, we found that the results derived here 
seem to generalize readily to configurations with a higher number of points. We found in 
all the cases studied, that if the viewed object is rotating with constant angular velocity 
around an axis perpendicular to the viewing axis, then the discrete Formulation yields 
a converging internal model, but that the rate of convergence depends on the angular 
displacements between consecutively viewed frames. 

In order to define the concept of rate of convergence, we first point out that recursive 
use of Eq. (3.21) shows that. f (A:) can be readily computed from <(()): 

k l 
r(A:)-(nB(j)k(0). (3.22) 

3 

and thus the convergence of <_{k) to depends only on this matrix multiplication, rather 
than on the perturbations themselves. Furthermore, if B(j) is cyclic (as it is for rotations 
in which the ratio between the angular displacement and 2tt is rational), with cycle m, 
convergence* depends only on the multiplications of the matrices over the cycle, which 
we denote: 

TO. 1 

B m = [] B(i). (3.23) 

)= o 

Interestingly, the behavior of B J m bears a strong similarity with geometric sequences. 
Let us define the; spectral radius of B m , p(B TO ), as the maximal modulus of the charac- 
teristic values of B m . Then it is possible 1 to show (see Appendix D) that: 

p?(B m ) = P (Bl). (3.24) 

From this result, one can understand the following important conclusion (see, for exam- 
ple, Varga, 19G2): 

lim Bj, ■-= *> ( >(B m ) < 1. (3.25) 

j -too 

In other words the necessary and sufficient condition for the internal model to be asymp- 
totically stable around the true solution, is that the spectral radius of the rotation matrix, 
B TO , is less than 1. 

Furthermore, from Eq. (3.24). one sees that the spectral radius corresponds to a 
"time constant" by which we can measure the velocity of convergence of the internal 
model. In fact, as p becomes closer to 1, more revolutions are necessary to obtain 
the same amount of convergence of the system, and the spectral radius of B^ r , declines 
exponentially with j. This exponential decline, however, holds only for the case of small 
perturbations, from which Eq. (3.22) was derived, and describes only a general trend, as 
in fine detail, the cyclic matrix B 7 „ is composed of a multiplication of partial matrices 
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that may impose 1 a "noise 11 in the decline. This noise can he seen in the simulations 
shown in section 3. 1. 

The dependence of the rate of convergence of the internal model on the angular 
displacement between viewed frames, d, and its consequences, are belter illustrated by 
some examples. Fig. 8 displays p as a function of d for the case in which the viewed 
object appears as an equilateral triangle, when viewed along the rotation axis. 




30° 60° 

Angular Displacement 



90° 



Figure 8. Graph of the spectral radius p as a function of d for the case in which the viewed 
object appears as an equilateral triangle, when viewed along the rotation axis. 



The first result of interest in this graph is that p < 1 for every angular displacement 
d such that 0° < d < 90°. This means that for almost every displacement used by the 
algorithm, the true solution (up to a depth reversal) is an asymptotically stable one. This 
result held for every configuration tried, provided that the rotation axis was perpendicular 
to the viewing axis. Thus, in spile of the fact that the true solution is asymptotically 
stable in the case of the discrete formulation, in the limit of arbitrarily closed frames, 
that is, in the continuous formulation, the true solution is not asymptotically stable. As 
can be seen in the figure, the; spectral radius tends to 1 as d — > 0. Thus, in the limit, 
the condition expressed in Eq. (3.25) for the convergence of B:J U to 0, is not fulfilled. 
This confirms our previous results regarding the lack of stability for the continuous 
formulation. 

In addition, from the limit of 1 at small displacements, p decreases with d, showing a 
steady improvement in the rate; of convergence of the internal model to the true solution 
with displacements up to 60°. This observation was also made- by Tillman (1984). The 
behavior of p with small displacements is examined more closely in the next section. In 
the rest of this section we discuss the affect of large angular displacements on the spectral 
radius. 

For large displacements, the rate of convergence deteriorates as can be seen in Fig. 
8, approaching the value of I for d — 90°. The spectral radius is 1 for displacements 
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of 90" because only two different views of I lie object are available under orthographic 
projection, which is too little information for recovering; structure from motion based on 
the assumption of rigidity. This result appears to be independent of other properties of 
the moving points. 

This behavior of/; with angular displacement, can be related to the convergence rate 
observed in computer simulations. This is illustrated in Fig. 9. We began with the 
set of three points, whose true coordinates (when projected onto the X — Z plane) lie 
around an equilateral triangle, as shown in Fig. 9a. The y coordinates of the points were 
nonzero but small. These points were then rotated around the vertical axis, and the 
error measure shown in Eq. (3.1) was computed, as a function of time. Fig. 9b shows 
an error plot for the case in which the angular displacement, between frames was 20", 
and the points were rotated for a total of 10 revolutions. The corresponding graph, with 
error plotted on a log scale, is shown in Fig. 9c. It can be seen that on a log scale, the 
long term convergence of the computed 3 D model is approximately linear. A measure 
of convergence rate can be given by the slope of the line that best iits this data, in a least 
squares sense. In Fig. 9d, the slope of this line is plotted as a function of the angular 
displacement between frames (the three point, configuration was rotated for a total of 10 
revolutions in each case). It can be seen that this graph is qualitatively similar to the 
graph of p shown in Fig. 8. 

3.2.5 Convergence Under Small Angular Displacements and Instability of 
the Continuous Formulation 

The question addressed in this section is, why does the rate of convergence fall with 
decreasing angular displacement between consecutive frames, eventually leading to in- 
stability for the case of the continuous formulation? One would certainly expect poorer 
corrections for the perturbations if smaller displacements are used, because less resolu- 
tion is available in the image data. That is, the signal to noise ratio is substantially 
reduced. Many more views, however, are seen in the case of small angular displace- 
ments, which might trade off with the deterioration in the data resolution. One must 
keep in mind that there are two types of noise that can influence the performance of the 
incremental rigidity scheme: (1) the noise from the image data, and (2) the difference 
between the model and the; true solution. In the computer simulations presented ear- 
lier, the image measurements were computed analytically, and hence had little error, so 
that the primary source of noise was the discrepancy between the computed model and 
true structure. This source of noise contributes to the deterioration of the solution with 
smaller angular displacements between frames. From the results in the previous section, 
we can conclude that in some sense, with decreasing displacement, the deterioration due 
to poorer corrections occurs at a faster rate than the improvement with the number of 
views. We explore this phenomenon in this section. 

Formally, we can express the above question as, why does the spectral radius p 
of the rotation matrix B„, increase; with decreasing angular displacement d, eventually 
tending to 1 when d — > 0? This is certainly a characteristic of the curve; of p vs. d 
curve if d is small enough (see Fig. 8). The matrix B m results from the multiplication 
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Figure 9. The results of computer simulations, (a) The true three point configuration, projected 
onto the X — Z piano, (b) The error in the computed 3D model, as a function of time, for 20° 
displacements between frames, (c) The error in the 3 D model, plotted on a log scale, (d) The 
rate of convergence of the solution, as a function of the angular displacement between frames. 
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of the matrices corresponding to partial displacements (Va\. (3.2'!)). One would like to 
define some quantity related to these partial displacement matrices, which expresses the 
deterioration of their informal ion content, and to study the connection of this quantity 
to the general spectral radius as a function of d. This quantity and its relationship to 
the spectral radius we're found by numerical experimentation. The intuition underlying 
these observations is what follows. 

(liven three point configurations and fixed angular displacements, an inspection 
of the partial displacement, matrices shows that they all have exactly the same two 
characteristic values, cj and c^-, which are real, and < cy < I and 1 < c-2- If these 
characteristic values belonged to a diagonal matrix, then the Hurt of such a matrix would 
be to change the perturbation vector ( -- (<i- (, 2) T into (rv l -<-2< 2)- hi general, however, 
the vector < can point in any arbitrary direct ion, and the effects of a diagonal matrix 
on the length of t would be different from direction to direction. Thus a meaningful 
quantity belonging to these matrices must be some average 1 of the effect of the matrix 
over all directions of the perturbation vector. An important example is the mean length 
squared of the vector resulting from the application of the matrix to the general unitary 
vectors, (cos^,sin.^). This quantity is the following: 

-L r^cWW + c2*m 9 ( 7 ))rf7 = ^4-- (3-26) 

2tx J 2 

The partial displacement matrices arc in general not diagonal, but' if d is small 
enough they can be closely approximated as orthogonally similar to the diagonal matrix 
having C\ and c-i as its elements. If we begin with an arbitrary perturbation vector, 
these matrices will rotate the vectors and transform the resultants as described above, 
and then repeat tins transformation until a cycle is completed. Thus for small enough 
displacements, one can expect the diagonal matrices to work on perturbation vectors of 
almost all directions, thus having an average effect as described in the last paragraph. 
Generally this mean effect always represents a decrease in the perturbation vector, which 
is applied 2tt/c/ times during a cycle of revolution of the object. Thus for small enough 
displacements one can expect a connection between ci, c%, p and d given by an equation 
of the following form: 



cl(d) + 4(d) 



lir/d 



p(d) = ^ni^Lj^Lj . (3.27) 

Eq. (3.27) closely matches the behavior of the spectral radius, for the equilateral 
triangle case (Fig. 8) for d < 20°. For other configurations, relationships very similar to 
this one were observed to apply. The importance of Eq. (3.27) is that it relates a quantity 
that can be computed in any partial displacement matrix, namely (c^(d) + c? i (d))/2 (the 
deterioration factor), to the spectral radius of the full rotation matrix, and it enables us 
to reduce the problem of why the continuous formulation is unstable to the question of 
how fast does the deterioration factor tend to 1 as d — > 0. In fact one can see from Eq. 
(3.27) that a necessary condition for liny > () p{d) — I is that,: 
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«( »o 2 ; 

This condition, however, is not sufficient. For example, if the deterioration factor 
approaches 1 linearly with d, i.e. (c'i(d) \ c\(d))/2 »i 1 ad for small d, then because of 
the power 2n/d in E<]. (3.27) we would have liny ,<> />(d) - c 2,r " < 1. We expert, then, 
that the deterioration factor approaches I faster than linear. This is indeed the case, as 
can be seen by expanding (c-\(d) I 1.2(d)) /'2 into the first few powers of its Taylor series. 
For any three point configuration, this expansion yields: 

iMLiW „ , . iU K (3.20) 

This quadratic approach of the deterioration to 1 is fast enough to ensure that liny ,o ~ 
1, as can be seen by substituting Eq. (3.29) into Eq. (3.27) and taking the limit. 

Therefore, we confirm that the instability of the continuous formulation is due to 
the fact that the deterioration in the data resolving power as a result of the small angular 
displacements between frames is faster (quadratic, with the inverse of the displacement) 
than the increase of information due to the increase in the number of frames (linear with 
the inverse of the displacement). This analysis provides theoretical support for the ob- 
servation that for the case of the discrete formulation of the incremental rigidity scheme, 
the rate of convergence of the internal model to the true solution diminishes as the an- 
gular displacements between frames decreases. In the limit of very small displacements, 
that is, in the case of the continuous formulation, this deterioration leads to a lack of 
asymptotic convergence of the model to the true structure. 

4. Perspective Formulations of the Incremental Rigidity Scheme 

In this section, we first present discrete and continuous formulations of the incremental 
rigidity scheme (hat use perspective projection of the scene onto the image plane. We 
then present a few examples of the results of computer simulations, in which the per- 
spective formulations were applied to sequences of points undergoing both pure rotation 
about a single axis in space and pure translation through space. The behavior of the 
perspective formulation of the incremental rigidity scheme is more complex than that 
of the orthographic formulation. We found that if the absolute position of an object 
in space is known throughout the motion of the object, then the perspective formula- 
tion performs well, similar to the orthographic formulation. If the absolute position of 
the object is unknown, the incremental rigidity scheme in general does not derive the 
correct 3 D structure. The results of computer simulations revealed a degradation in 
performance with smaller angular and spatial displacements between frames, but this 
degradation was somewhat more severe than in the orthographic formulation. In the 
limit of the continuous formulation, the solution is no longer stable. 

In both the discrete and continuous formulations, we assume that the positive z axis 
points in the direction of the optical axis, with the image plane at z — 1, as shown in Fig. 
10. For notational convenience, we let {xi(l-),ffi(t),Zi(t)) represent the true coordinates 
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of point: i in space and («,:(/•), '",:(')) denote the (rue projected image coordinates of point 
i. The image coordinates arc given by: 






(4.1) 



The coordinates of the points in the computed model are denoted by (xj(t),yi(t), Zj(t)). 
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Figure 10. The imaging geometry for perspective projection. The positive z axis points in the 
direction of the optical axis, with the image plane at z = 1. 



4.1 The Discrete Formulation 

In the discrete case, we compute the depth values Zi(t') that minimize the measure of 
rigidity, D,i{t,t') given by the following (for notational simplicity, (a:;,//,-,,?,-) refer to the 
coordinates of the points at time t in the computed model, while (^!-, ?/',£;) refer to the 
coordinates of the points at a later time t ! ): 

= E [(( x >- - x ^) 2 + (»» - ?/;) 2 + te ~ *>) 2 ) * (4-2) 

3 

((^■-*;) a + w-»;) v K^-4) a )* 



*>j 
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This measure is (he same as (hat used in (lie orthographic formulation (i.e., it is derived 
from Eq. (2.1) with the I'jAl) factor in the donomonator omitted, as shown in Eq. (Al) 
of Appendix A). 

As before, we assume that we have a current model (x;(/,), ?/,;(/.),£,;(£)). Rather than 
assuming' that. Xi{t') and y,(l') are known, however, we assume thai, only «,-(/.') and t>i(t') 
are known (they are derived from the true s]>are coordinates at time /.'). The coordinates 
of the points in space at time /.' for the computed model are then given implicitly by: 

Xi (t') = Uiit'W) 

15 y substituting Eqs. (4.3) info Eq. (4.2), D,i(l,l') can be expressed in terms of the 
coordinates of the points as follows (('«[-, v\) denote the image coordinates at. time £'): 






(4.4) 



Aft(T the values •?;(£') that minimize Dj(t,t') are computed, Eqs. (4.3) are used to derive 
the space coordinates x.-(t') and yi{t'). The new 3D coordinates, (z; (£'),;/;(£'), ^i(i')), 
then become the current 3 D model, a new frame in the sequence is registered, and the 
process repeats itself. 

In the case; of perspective projection, 3 D structure can be recovered from relative 
motion only up to a multiplicative scale factor. This is clear by inspection of Eqs. (4.1); 
a constant scaling of the space coordinates {x-i{t) , yi(t) , Zi{t)) does not change the image 
coordinates (v,i(t),Vj(t)). In most of the computer simulations, we assumed that the 
overall scale factor is known; this is equivalent to assuming that the absolute position of 
the object in space is known. In the simulations, the absolute coordinates of one of the 
n points in motion was supplied to the algorithm and the coordinates of the remaining 
n — 1 points were computed relative to the known point (similar to the simulations of 
the orthographic formulation). If the absolute position of the object in space is not 
known, then the scale of the computed 3 D model depends on the choice of the initial 
configuration. For example, suppose that we assume that the set of moving points 
initially lies on a plane (hat is parallel to the image plane. The scale of the computed 
3 D model will then depend on whore this initial plane is placed in depth. That is, if 
the plane; is positioned twice as far from the image plane, the computed 3 D model will 
essentially be twice as large. We will briefly mention the results of computer simulations 
in which the absolute position of the object in space is unknown. In these simulations, 
when we began with a flat; initial configuration, we usually placed the initial z coordinates 
of the points at a depth that was approximately the mean of tin 1 true z coordinates of 
the points in space. Further details of the implementation of the discrete formulation 
can be found in Appendix A. 
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4.2 The Continuous Formulation 

In the continuous case, wo compute the z components of velocity that minimize the 
measure of rigidity, D,-{1), given by the following (the coordinates and velocities are all 
measured at time /.): 



v^ \( x i - x jW r -i - *j) + hn - Vj){.vi - in) •+ ( g » - z j)i i j r3')l. a 

(x.i - x.j) 2 f ('//,• - ;</y) 2 4- (z; - ;?y) 2 



Dc(i ) ,, £ I^IL 1^2 zr ^^_^^i^ T _^^^ T ^l^i^^i_. (4 . 5) 



This measure is the same as that used in the orthographic case (see Eq. (2.7)). It is 
assumed that a current model of the coordinates of the points in space is given, along with 
the known coordinates and velocities of the points in the image plane, ».,;(/.). w,(/.), «;(/•), »i(t). 
The velocities of the points in space for the computed model are then given implicitly 
by: 

Xi(t) =Ui(t)Zi(t) + Ui[t)Zi(t) 

ili{t) = Vi{t)zi{t) + vi{t)zi{t). (4.6) 

If we let Xij — xi - Xj, y.;.j = yi - j/y, 2,-y = Zi - 2y, and % = £; - iy, then Z? c (i) is 
ex{>ressed in terms of the coordinates of the points as follows: 

n a\ - V [ x *j( M » ir ' + iOHi ~ u Ai ~ il -j z j) + VijbHZi + "»-g t - - vjzj - %fy) + %%] 
,:,y ^'i H % + *tf 

(4.7) 
In other respects, the continuous formulation is similar to the discrete formulation. A 
model of the structure of the moving object is built up by continually taking into account 
new velocity information over an extended time period. After the Zi{i) are computed 
by minimizing Eq. (4.7), Eqs. (4.G) are used to derive X{(t) and y%(t), which are then 
used to compute the new 3D model. Again, because perspective projection is used, the 
structure of the object can only be computed up to a multiplicative scale factor. Further 
details of the implementation of the continuous formulation are presented in Appendix 
A. 

4.3 Computer Simulations 

In this section, we present some results of the computer simulations of the discrete 
and continuous perspective formulations of the incremental rigidity scheme, for a small 
number of examples. We first consider the case in which a set of discrete points is 
translated through space, and then consider the rotation of a set of points about a 
vertical axis. Through simulations of the discrete formulation, we examine the rate of 
convergence of the algorithm and the quality of the solution that it yields, as a function 
of the size of the angular and spatial displacements between frames. We then observe 
the limiting behavior of the continuous formulation for the case of rotation. We also note 
the difference in performance of the perspective formulations when the absolute position 
of the moving points in space is either known or unknown. 
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K should be stressed that it) lite case of orthographic projection, it- is not possible 
to recover 3 1) structure for objects undergoing pure translation, because (lie relative 
positions of projected points in the image do not change as objects translate. The use of 
perspective projection has the advantage of allowing the recovery of structure for trans- 
lating objects. In the simulations presented here, objects were oscillated back and forth 
in the direction parallel to the image plane. If the objects were allowed to translate 
hi one direction for a large extent, the projected image of the points would eventually 
become so small that the recovery of structure would become- difficult due to a loss of 
numerical accuracy in the measurement of the changing positions of the points. Oscil- 
lating the points allows us to analyze their structure over an extended time period while 
maintaining a relatively large image. 

In the first set- of examples, the input consisted of a set of six points in space. The 
coordinates of the points were chosen randomly and the coordinates of oik 1 of the points 
was known to the algorithm. The projections of the points onto the X - Y and X — Z 
planes are shown in Figs, 11a and lib, respectively. This set of points was oscillated back 
and forth in a direction parallel to the image plane. The X — Z projection (birds' eye 
view) of the rightmost and leftmost position of the points in space are shown in Figs, lie 
and lid, respectively. The x coordinates of the points in Figs, lie and lid are shifted 
by GO units from the initial positions shown in Fig. lib. In the simulations, the initial 
points shown in Fig. lib were first translated to the right in constant discrete steps and 
then translated to the left with the same discrete steps. A 3 D model was built up over 
several oscillations of the points. The initial model for the structure of the points was 
flat, with the z coordinates placed near the mean of the true z coordinates of the points. 
For the particular sot of points used in these example's, Zi(t) = 80, i — l,...,n, in the 
initial 3D model. Once the initial Z{(t) are specified, the initial Xi(t) and yi(t) can be 
determined from the image coordinates Ui(t) and v t {t). The one point whose position is 
known througlunit the oscillation of the points is indicated by an arrow in Fig. 12a. 

Fig. 12 shows a birds' eye view of the computed 3-D model after 1 (left column), 
4 (middle column) and 10 (right column) complete oscillations of the points, using three 
different sizes of spatial displacements between frames. The true 3D structure of the 
points is shown again in Fig. 12a. In the case of Fig. 12b, the size of the spatial 
displacements between frames was GO, so that each complete oscillation consisted of four 
frames (i.e., the points were displaced to the right in one step, returned to the central 
position, displaced to the left in one step and then returned again to the central position). 
In the case of Fig. 12c, the points were translated in steps of 30, so that a complete 
oscillation of the points consisted of 8 frames. Finally, in the case of Fig. 12d, the points 
were translated in steps of 5, so that each full oscillation consisted of 48 frames. In Fig. 
12e, we show plots of the error in the computed models as a function of time, for 10 
complete oscillations of the points. The error measure xxsed here is the same as that used 
in the orthographic case, shown in Eq. (3.2). The error plots for spatial displacements 
of GO, 30 and 5 are shown superimposed, with the displacements indicated above each 
curve. From the results shown in Fig. 12, it can first be seen that the algorithm does 
essentially converge to the true structure of the points. Second, the rate of convergence 
of the computed model to the true solution decreases with smaller spatial displacements. 
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Figure 11. A set of six points with random coordinates is shown projected onto (a) the X — Y 
plane and (b) the X — Z plane, (c) The rightmost position of the points during tlieir oscillations 
hack and forth in the direction parallel to the image plane, (d) The leftmost position of the 
points. 



For example, after a single oscillation (leftmost, column), the model shown in Fig. 12b 
(displacements of GO) is clearly closer to the true structure of the points shown in Fig. 12a 
than the model shown in Fig. 12d (displacements of 5). The three error plots shown in 
Fig. 12e also have periodic variations superimposed on a steady decline in error, with the 
amplitude of the variations increasing with smaller spatial displacements. Tims, in the 
case of tin; perspective formulation, it appears that the performance of the incremental 
rigidity scheme degrades with smaller spatial displacements when objects undergo pure 
translation through space. This result is analogous to the observation that in the case of 
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orthographic projection, the performance of the scheme degrades with smaller angular 
displacements when objects undergo rotation. 
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Figure J 2. (e) Graphs of the error in the distances between points in the computed 3 D model, 
as a function of time. The points underwent JO complete oscillations, in steps of GO, 30 and 5. 



Fig. 13 illustrates the behavior of the perspective formulation, for the case of rotation 
of the points about a central vertical axis. The set of six points shown in Fig. 1.1 was 
rotated about a seventh point that was placed at the position: 



(x i {t),y i {t),z i {t)) = (0,0,80). 

The axis of rotation passed through this point and was perpendicular to the X- ~ Z plane. 
The position of this central point was known to the algorithm throughout the motion of 
the points and the positions of the remaining six points were computed relative to this 
central point. As before, the initial configuration was flat, with Z{(t) = 80, i — l,...,n. 

The set of 7 points was rotated in angular steps of 40° and 5°. Figs. 13b and 13c 
show the computed 3 D model after total rotations of 40°, 80°, 120° and 3G0° (indicated 
on the left), for angular displacements between frames of 40" and 5°, respectively. The 
true positions of the points are shown in Fig. 13a. In Fig. 13d, the graphs of the error 
in the computed models are shown as a function of total angular rotation, for three 
full revolutions of the points. It can be seen that the rate of convergence and quality 
of the 3D structure decreases with the smaller angular displacement between frames. 
Although similar to the case of orthographic projection, the degradation in performance 
with smaller angles was somewhat more severe in the case of perspective projection. 

Fig. 14 illustrates the behavior of the continuous formulation of the incremental 
rigidity scheme, for the case of perspective projection. The same set of points shown 
in Fig. 11 was rotated about the vertical axis and the 3D model was computed at 
infinitesimally closely spaced frames, using the instantaneous velocities projected onto 
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Figure 12. (a) Projection onto the X — Z plane of the true 3D structure of the set of points, 
(b) The computed 3D model sift or 1 (left column), 4 (middle column) and 10 (right column) 
complete! oscillations of the points, for a spatial displacement between frames of size GO. (c) 
and (d) The same three computed 3-D models, for spatial displacements of size 30 and 5, 
respectively. 
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the image. As in the previous rotation examples, a seventh point was placed at. the renter 
of rotation of the points, and (lie position of this central point was known to the algorithm 
throughout the rotation of the points. In Fig. Ha, we compare the true positions of 
the points (left) with the best solution (right) obtained over four lull revolutions of the 
points. Although the model is (mite close to the true structure at, this position of the 
points, the solution oscillates significantly over an extended time period, similar to the 
orthographic case. A graph of the error in the computed model over the four revolutions 
is shown in Fig. 14b. The arrow marks the point at which the solution shown in Fig. 14a 
was obtained. There is an initial fast convergence toward the true structure of the points, 
but the algorithm then oscillates with high amplitude away from the true solution. We 
conclude that similar to the orthographic case, the direct use of velocity information for 
the recovery of structure by the incremental rigidity scheme does not lead to a robust 
solution. 

We should finally note that we also explored some examples in which no information 
about the absolute position of the points in space is known. In this case, we began with 
an initial configuration that was Hat, and otherwise provided no further constraint on 
the positions of the points. The points were either oscillated back and forth in the 
direction parallel to the image plane, or rotated about a vertical axis. We found that 
the incremental rigidity scheme would sometimes derive the correct solution under these 
conditions, but in general, the algorithm did not derive the correct structure, although the 
computed solution was essentially rigid. The computed solution was not just a scaled 
or rotated version of the true structure, but actually a different one altogether. This 
suggests that some additional constraint may be required to derive the correct solution. 

To summarize, the perspective formulation of the incremental rigidity scheme ap- 
pears to perform well when the absolute position of the object in space is known. This 
formulation then allows the recovery of structure for objects undergoing pure translation 
through space, as well as rotation. The results of computer simulations revealed a degra- 
dation in performance with smaller angular and spatial displacements between frames. 
This degradation appeared to be more severe than in the case of orthographic projection, 
when the points were; rotated about a vertical axis. When the points were translated, 
there was a more gradual decrease in the rate of convergence and quality of the computed 
3-D model with smaller spatial displacements. In the limit of the continuous formulation, 
the solution was no longer stable. 



5. Summary and Conclusions 

In this paper, we studied and generalized the incremental rigidity scheme for the recovery 
of structure from motion. This algorithm, as first proposed by Tillman (1984), assumed 
the visual input to consist of a sequence of frames, each containing a finite number of 
points. The scheme maintained an internal model of the structure of the viewed object, 
which was updated from frame to frame, so as to be consistent with the changing image 
and to be as rigid as possible. This internal model was shown to correctly converge to 
the true structure of the object, for both rigid and nonrigid motions. 



39 



40 



Z i 

« " » 

X 



Z" 



z 



30 



Z" 



Z ■' 



Z" 



120 



Z" 

• e 

* « 

X 



z- 



360 



Z 



Z- 



Figure 13. (a) The true positions of the points after rotations by 40°, 80°, 120° and 3G0°. 
(b) The computed 3 D model after rotation of the points by 40°, 80°, 120° and 300°, for an 
angular displacement between frames of 40°. (c) The same computed 3 D models, for an angular 
displacement of size 5°. ((d) on next page.) 
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Figure 13. (<1) Graphs of the error in the distances between points in the computed 3 D model, 
as a function of time, for three full revolutions of the points. 



Ill the present paper, we focused on an observation made by Ullinan, that the rate 
of convergence of the internal model to the true structure increased with the amount 
of change between consecutive frames. A major part of our analysis focused on objects 
that rotated with constant angular velocity around an axis perpendicular to the viewing 
axis, and whose image was formed through an orthographic projection, at equally spaced 
angular displacements. The; complete dependence of the convergence rate on the size of 
the angular change was obtained for a large variety of objects. 

Tillman's observation was confirmed for small displacements. Indeed, the conver- 
gence rate first increased towards a maximum and then decreased, as a function of the 
angular displacements between the frames. The deterioration of the algorithm for high 
displacements can be understood if we note that the number of frames available for anal- 
ysis in a given total number of revolutions of the object, or in other words, the amount 
of information provided to the scheme, is reduced in inverse proportion to the angular 
changes between the frames. 

The deterioration of the algorithm for small displacements, however, has a more 
complex and surprizing basis. In this case, although the number of frames used in a 
given amount of rotation is high, the spatial resolution between the points is low (see 
discussion on the beggining of section 3.2.5). Our analysis indicates that the latter effect 
dominates for very small angular jumps between frames, reducing the convergence rate. 
Furthermore, we found that this reduction is such that in the limit of infinitesimally 
small displacements, where only velocity information is used, the internal model does 
not even have a full convergence to the true structure of the viewed object, but only a 
rough approach to the correct solution. 

In analogy to linear algebra, one can say that the transformations from frame to 
frame of the internal model in the incremental rigidity scheme become ill conditioned as 
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Figure 14. (a) The best solution (right) obtained over four revolutions of the points is compared 
with the true positions (left), for the case of the continuous formulation, (b) The error in the 
internal 3D distances between points in the model, as a function of time. The arrow in (b) 
indicates the instant at which the solution in (a) was obtained. 



the angular displacements decrease, i.e. they become sensitive to noise. This property 
is not unique to the present transformations, but occur in other important situations 
such as numerical differentiation. The noise sensitivity of a differentiation of order m 
increases as n m when the number of points n in a given interval is increased (Strang, 
197G). 

As mentioned before, the orthographic, projection is in general not physically valid. 
It was, used in the present paper, because it allows a deeper analysis of the phenomena 
studied. Tins analysis is further validated because these results are not limited to the 
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case of orthographic projection, or to the case in which the object, rotates. Under certain 
conditions, the incremental rigidity scheme can be generalized successfully to formula- 
tions that, use perspective projection, and that in this case, the structure of the viewed 
object is recoverable both for rotation and translation. Also, similar to the orthographic 
case, the convergence rate of the internal model to the true st ructure decreased, if smaller 
spatial changes between consecutive frames were used. 

Inspection of previous results by other researchers suggests that the use of informa- 
tion that is local in space and time is insufficient to allow a robust- recovery of structure 
from motion. We distinguish between two types of local information that, are relevant to 
the problem: (I) spatial locality, which refers to the use of a small number of points or 
features of the viewed object,, and (2) temporal locality, which is the use of a small num- 
ber of views of the object. Early algorithms for recovering structure from motion based 
this recovery on a limited number of points .and views, and were able to recover the struc- 
ture of the viewed object, when it was rigid. When an object deformed, some algorithms 
could determine its nonrigidity through inconsistencies of the underlying equations. It 
followed that these schemes are not robust against noise, as noise has similar effects as 
nonrigidity. 

It. can be shown that motion information that, is extended in space, but temporally 
local, can be used to provide stable recovery of structure from motion (Brass and Horn, 
1983; Negahdaripour and Horn, 1985). Perceptual evidence, however, suggests that 
spatial extension is not necessary for the solution of the structure from motion problem 
by the human visual system (Borjesson and von Hofsten, 1973; Johansson, 1975; Petersik, 
1980). 

On the contrary, Ullman's incremental rigidity scheme, which was applied locally in 
space, both in its original formulation and in the present paper, uses an internal model 
that, is updated constantly, in order to extend motion information over time, and converge 
to the correct solution. Some perceptual studies suggest that this temporal extension 
may be used by the human visual system (Wallach and O'Conncll, 1953; White and 
Mucser, 19C0; Green, 1961). 

We showed in the present paper, however, that, an algorithm that, overcomes tempo- 
ral locality does not. necessarily provide a robust solution to the structure from - motion 
problem. It seems that in this case, a necessary condition for a stable solution is the 
use of views of the object that are significantly disparate. In the limit, the use of pure 
velocity information allows only a rough solution to the structure from motion problem, 
which is less stable over an extended time period. A somewhat similar conclusion was 
derived by Ullman (1983), who argued that in a temporally local scheme, the recovery 
of structure from the instantaneous velocity field is impossible under orthographic pro- 
jection, and that for perspective projection, the recovery is unstable. We suggest that 
the problems with using local velocity information are too large to be overcome by an 
extension of this information over time. In other words, the inaccuracies introduced by 
a velocity based scheme cannot, be corrected by the use of multiple views of the object. 

Our results therefore hint, that even in a temporally extended algorithm, a memory 
of sufficiently distinct past views or internal models of the viewed object is necessary for 
a robust recovery of its structure. We point, out, that the incremental rigidity scheme, in 



43 

its original formulation and in the present paper, used memory of only one past, internal 
model. Although in principle, a memory of many past models could be used by the 
algorithm, it is remarkable that such a minimal amount of past information is sufficient 
to provide 1 a robust recovery of structure. Still, it may be \]\v case that, the use of multiple 
past views or internal models for establishing a new model of 3 D structure could increase 
the rate of convergence of the internal model. 

hi further support of the above* conclusions, we note that in order to recover 3 D 
structure from measured image velocities, it is essential that these velocities at least ap- 
proximate the true projected velocities of moving elements in space. In general, however, 
it is difficult to derive real projected velocities from the optical flow pattern on the eye 
or camera. Factors such as changing illumination, spccularities, shadows, and rotation 
of an object, surface relative to a light source, can create patterns of optical flow that do 
not correspond to the real movement of features on a physical surface. The difficulty of 
computing correct, projected velocities provides an additional motivation for not basing 
the recovery of structure from motion directly on velocity information. 

Our analysis raises a number of issues regarding the recovery of structure from 
motion in the human visual system. First, it is not clear whether the visual system 
achieves a stable solution to the structure from motion problem, or a rough solution 
such as that provided by a velocity based scheme. The more robust solution may not be 
essential if we consider that other perceptual cues, such as binocular stereo, may help 
to improve the quality of the 3D solution obtainable; by the structure from motion 
process. Further psycophysical experiments are required to examine this question. 

A second point of interest is that our analysis suggests that if the visual system 
incorporates a robust solution to the structure -from motion problem, it must be able to 
match corresponding points in very disparate views of moving objects. The displacements 
between corresponding points may be larger than the spatial limits proposed for the 
short-range perceptual process found for 2 D motion (Braddick, 1973, 1974, 1980; Pantle 
and Picciano, 1976; Petersik, Hicks and Pantle, 1978; Petersik and Pantle, 1979). A 
similar conclusion, based on different grounds, was formulated by Petersik (1980). He 
explored the sensations elicited by stroboscopic simulations of a rotating transparent 
sphere filled with randomly positioned dots. By manipulating both spatial and temporal 
variables in the simulation, he conchided that, corresponding elements in consecutive 
frames can be matched over spatial and temporal distances that exceed the empirically 
determined limits of the 2 D short range process. Even in this experiment, however, 
points that reach the periphery of the sphere have small displacements from frame; to 
frame that may fall spatially inside the short, range process. It remains, therefore, to be 
established that; these points do not proviele all the information used by the subjects to 
sense the continuous rotation and internal volume of the sphere. 

If the human visual system is able to match very distant cenrcsp ondmg points from 
disparate views of a moving object, the question e>f how we solve this correspondence 
problem becomes important (for a review of the motion correspondence problem, se;e 
TUliium (1981)). This correspondence could be established tnthor through the repeated 
use of a short range matching process, or through the use of an explicit lemg range 
matching process. In the first e'ase', the short range process could provide an essentially 
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continuous (racking of features in (lie changing image, and (lie positions of features could 
then he sampled by a longer range tracking process that feeds disparate positions into 
the structure from motion process (a similar idea has been suggested by Ullman (1981)). 
Iu the latter case, a correspondence would be established directly from disparate views 
of moving objects, without the aid of the intermediate short range process. Petersik's 
(1980) results may suggest that it is possible to solve this long range correspondence 
problem for recovering structure. This correspondence problem could be solved in con- 
junction with the structure from motion process; that is, a correspondence could be 
chosen that subsequently gives rise to the most rigid 3 D interpretation (such a match- 
ing process requires a global decision procedure, and may be well suited to solution by 
parallel schemes such as Tlopfield's "neuronal nets" ( J 98d ) ) . Given the difficulty of solv- 
ing this long range correspondence problem in general, however, it seems unlikely that 
short range measurements would not be used when they are available. 
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Appendices 

Appendix A 

In t.li i s appendix, we present some of the details of the implementation of the incremental 
rigidity scheme used to derive the simulation results presented in sections 3 and 4. 

Orthographic Projection: The Discrete Formulation 

Ulhnan formulated the incremental rigidity scheme as the computation of depth values 
Zi(t') that minimize the measure of rigidity given hy D{t,t'), as shown in Eq. (2.2). 
The analyses presented in this paper mainly consider rigid objects in motion, which are 
compact in the sense that the internal distances between pairs of points do not differ 
much from one another, hi this case, the additional distance factor in the denomonator 
of the measure used by Ulhnan has little influence on the resulting solution obtained 
by the algorithm. We therefore removed this factor in most of our simulations, and 
minimized instead the following expression (the distances Uj{t) and hj{t') are expressed 
in terms of the coordinates of the points, and (x,i,yi,Zi) refer to the model coordinates 
at time t, while [x'^y'^z^) refer to the model coordinates at time t'): 

Da(t, t') = J2 [((*< - *;) 2 + fa " %') 2 + (* - *i) 2 Y 

The Zi{t') that minimize D,i(t,t') satisfy a system of nonlinear equations given by: 

£|L =0 , , = !,...,„-,. { A2) 

Rather than solving this system of equations explicitly, we solved them implicitly through 
the use of a steepest descent minimization algorithm for Eq. (Al), based on the gra- 
dient of D,i(t,t') (Luenberger, 1973). The components of the gradient arc given by the 
following: 

In the case of orthographic projection, only relative depths can be recovered. In the 
implementation, we placed the initial point (xu{t),yu(t),Z\)(t)) at the origin of the co- 
ordinate system; that is, (xu(t).yu(i), z${i)) — (0,0,0). At every instant, the remaining 
n — 1 points were then given relative to this initial point. Unless otherwise stated, the 
initial configuration of points was flat, that is, Z{({)) = 0, for i,...,n — 1. From this 
configuration, the algorithm can move toward two equally likely configurations, one be- 
ing the mirror image reflection of the other about tin 1 image plane'. In other words, this 
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configuration is located at. a saddlepoint in the solution space for the problem (the solu- 
tion space* is an n dimensional space in which a value for D,i(t, t') is assigned to every 
possible combination of depth values z(t')), at which the gradient of D,i(t,t') is zero. To 
make this gradient nonzero, the initial solution was perturbed slightly, thus causing the 
algorithm to move in a particular direction. 

The steepest descent- minimization met hod from nonlinear programming (see, for ex- 
ample, Luenberger, 1973) was used to compute the Zi[t') for each new frame*. The current 
depth values Z{(t) were used as the initial solution for Zi(t') to begin the minimization 
for the next frame. The objective function is given by D,t(t,t'), and its gradient by Eqs. 
(A3). Golden section search was used to perform the one dimensional minimization 
within the steepest descent method. 

Orthographic Projection: The Continuous Formulation 

In section 2, we presented a continuous formulation, which requires the computation of 
z components of velocity, i,-(f) that minimize the measure of rigidity given by D c (t), as 
expressed in Eq. (2.3). In the computer simulations, we used the particular measure of 
overall deviation from rigidity given by: 

n (t\ - V^ [( ^ - x 3 -){xi - ±j) + (yj - yj){in - j)j) + [zj - Zj){zi - Zj)} 2 
aJ ~ 2 -' (s.--*i) 2 + (l/i-yy) 2 + (2.--*y) 2 ' ( } 

The Zi[t) that minimize Da{t) satisfy a system of linear equations given by: 

dDr 0, t = l,...,n-l. (A r >) 



dZi(t) 

These equations are of the form: 

dDc _ „ y- (xj - Xj)(xj - ij) + [yj - yj){iii - yj) + (g,- - Zj){zj - zj) _ 

dzi l 2-< ~ te-xyj' + fa-yy)' + (*-*,-)» {Zi ~ Zj) U - 

(AG) 
In the case of orthographic projection, only relative z components of velocity can be re- 
covered. As in the discrete case, we again placed the initial point {xa(t) ./y (t) , zo(t)) at the 
origin of the coordinate system throughout the entire motion, so that (xo(t),y^(t), Zo(t)) — 
{xQ(t),yo(t),zo(t)) = (0,0,0). At every instant the remaining n — 1 z components of ve- 
locity were then given relative to this initial point. Unless otherwise stated, we again 
began with a flat initial configuration in which .c'i(O) = 0, for l,...,n — 1, which was 
perturbed slightly so that the gradient of D c (t) is initially nonzero. 

At each moment, Z{(t) were obtained by solving a system of linear equations, for 
which we used the simple Gauss -Seidel relaxation method. The initial condition for the 
relaxation was usually the set of velocities computed in the previous iteration and i; = 
for the first iteration. To integrate motion information over an extended time period, we 
then, made use of the following approximations using i,-(<), iji{t), and i,-(<): 
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Xi(t) f Xi(t)U = xi[t \-fit) 

Vi{i) I !)i(W ---- Vi{t + fit) (AG) 

Zi(t) + ki(t)H = Zi(t + U). 

Once the Zi(t) wore computed, by minimizing D c (l), then (lie £,;(/ + fit) could bo derived 
from the current, model 3,-(i). The new model was then taken as M(£ f M) = (x{(t + 
fit),yi(l -f fit),Zi(t + <*>/,)), and the process was continued. The time interval, #£, typically 
corresponded to an angular displacement between frames on the order of 0.1 degrees. We 
also experimented with angular displacements up to two orders of magnitude smaller, 
and found the qualitative behavior of the algorithm to be the same. 

Perspective Projection: The Discrete Formulation 

In the case of perspective projection, the x and y coordinates at time /,' can be expressed 
in terms of the known image coordinates at time t' and the depth values to be computed, 
Zi(t'), as follows: 

Xi(t') = Ui{t')zi(t') 

The measure of rigidity to be minimized is given by substituting Eqs. (A7) into Eq. 
(Al): 

i 

2 



D d (t, t') = £ ((*.- - x i) 2 + {Vi ~ ?/y) 2 + (* - z i?) ■ 

l '> (A8) 

- {{u[z\-u)z)? + {v\z[-v>ztf + {z[--z> j )*y]\ 

The Zi(t') that minimize D c i(t,t') satisfy a system of nonlinear equations given by: 

° }Jd 0, • = l,...,n-l. (A9) 



dZi{?) 

In the computer simulations, these equations were solved implicitly through the use of a 
steepest descent minimization algorithm for Eq. (A8), based on the gradient of Dd(t,t'). 
The components' of the gradient are given by the following: 

The computed £{(£') are used to derive the new Xi(t') and yi(t') using Eqs. (A7), which 
then provide the new computed 3- D model. This new model also serves as an initial 
solution to begin the minimization for the next frame. 
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Perspective Projection: The Continuous Formulation 

In the continuous case, the x and y components of velocity of the points in space can be 
expressed in terms of the known image coordinates and velocities as follows: 

Xi{t) - Ui{t)ii(t) -\ Ui(t)Zi{t) 

Vi(t)^vi(l)zi(l) + Vi(t)zi(l). (All) 

If we let Xij = x,; — Xj, yij = ■}/{,— y 3 -, z,;j = Z{ - zj, and «,-y = i t - - iy, then D K (t) is 
expressed in terms of the coordinates of the points as follows: 



c ^ ' ■<£_/ .,-2 i ,,,2 i ,2 



• l 2 



»J 



I J -'tj I J 



(A 12) 



The ii(£) that minimize D,.(t) satisfy a system of linear equations given by: 



^ = » = l,...,n-l. (413) 

OZi 



These equations are of the form: 






dz 

3 



where the 6,-y arc given by: 



_ Xij(uiZi + iiiZi - UjZj - itjZj) f Vij{viZi + v-iZi - Vj ij - vj Zj) + Zijkij 

X ij ^ Vij ^ Z ij 

At each moment, the Z{(t) were obtained by solving the system of linear equations given 
by Eqs. (A14), for which we used the simple Gauss Seidel relaxation method. After the 
Zi(t) were known, the Eqs. (All) were used to derive the x and y components of velocity 
in space, i,-(i) and yi(t). To integrate motion information over an extended time period, 
we then made use of the approximations given in Eqs. (A6) to compute a new model 
{xi{t + 6t),yi(t + 6t),Zi(t + 6t)). 

Appendix B 

In this appendix, we explain how Eq. (3.10) was derived. 

Eq 3.8 depends on the image coordinates (data) and on depth coordinates of the 
points z.;(i) and their velocities Z{(t) (computed values). These depth variables are in 
general displaced from the true; motion values. This departure can be incorporated in 
Eq. (3.8) as follows: 

3D 

— 7 {t,zx +f.i,...,£ n _i +£„ -i,zi +ei,....i n -i -|-e„ i) = 0. (£1) 

OZi 
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If the c and < are small enough, as it is assumed in the stability analysis of section 
3.2.2, then the first elements of the Taylor expansion of Eq. (Bl) yields the following 
approximation: 

dD. „ . .. • 

=-7cr-_(t, *i, ..., z n . i,zi, ..., z n ^i) 

(j Z{ 

XT dD u - 

+ 2^^^~{^^-^ n .-uz l ,...,z n . l )c j 

3-1 

V± dD . . . ; , .. 

j - 1 l 3 
But the true rigid solution certainly minimizes Eq. (3.7) as it sets it to 0, thus being a 
solution of Eq. (3.8). It follows that the first term on the righthand side of Eq. (B2) is 
0, from which Eq. (3.10) is concluded. 

Appendix C 

In this appendix we indicate the method by which equations of the type 3.14 and 3.15 
were derived. 

Equation 3.14 is derived by taking the necessary partial derivatives of D (Eq. (3.7)) 
by Z{ and Z{ and combining them as indicated in Eq. (3.11). These derivatives cire 
evaluated at the true solution itself, that is at Z{(t) and Zi{t), and because we only 
considered rigid motions in this paper, we could use the following relationship to simplify 
the results: 

atj + (z{ - Zj){zi - Zj) = 0. (CI) 

Eq. (CI) is derived by setting the distances between points i and j as constant and 
taking the temporal derivative of this distance. 

The conclusion of Eq. (3.15) is derived by noting that, in a rigid rotation around 
an axis perpendicular to the viewing axis, with constant angular velocity w, zi can be 
written as: 

Zi(t) = eicos(iot + <f>i). (C2) 

Then, by direct substitution of Eq. (C2) in Eq. (3.14) and its integration as indicated 
in Eq. (3.15), we derive the stated result. 

The calculations explained in this appendix, though straightforward, are cumber- 
some, and were done by using Macsyrna, a computer system for performing algebraic 
manipulation. 

Appendix D 

In this appendix we show that if C is an n x n matrix over an algebraically closed field 
then: 
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fP(C) = p(C'), (Dl) 

where p is the spectral radius. 
Let, P be such that: 

P 1 CP = J, (D2) 

wliere J is in Jordan canonical form. It follows that: 

& = PPP- 1 . (J53) 

Similar matrices have similar characteristic values and thus similar spectral radii, thus: 

p(C) = p(3), (DA) 

and 

p(&) = p(J>'). (D5) 

Direct inspection of the matrices J 3 show that they have the same type of triangulariza- 
tion as J (that is, upper or lower triangularization), with the elements in the diagonal 
being the j - th power of those of J. Thus it follows that: 

p(3>')=p?(3), (DG) 

which together with Eqs. (C4) and (C5) imply the result stated in Eq. (Dl). 



